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SUMMARY 

The advantages inherent in the Boundary Element Method (BEM) for potential flows are 
exploited to solve viscous How problems. The trick is the introduction of a so-called Dual Re- 
ciprocal technique in which the convective terms are represented by a global function whose 
unknown coefficients are determined by collocation. The approach, which is necessarily it- 
erative. converts the governing partial differential equations (PDEs) into integral equations 
via the distribution of fictitious sources or dipoles of unknown strength on the boundary. 
These integral equations consist of two parts. The first is a boundary integral term, whose 
kernel is the unknown strength of the fictitious sources and the fundamental solution of a 
convection- free flow problem. The second part is a domain integral term whose kernel is the 
convective portion of the governing PDEs. The domain integration can be transformed to 
the boundary by using the Dual Reciprocal (DU) concept. The resulting formulation is a 
pure boundary integral computational process. 

INTRODUCTION 

The major advantage the BEM approach enjoys over other techniques is the confinement 
of the computation to the boundary. The result is the reduction in the effective dimension 
of the problem. The efficiency with which linear problems in rontinua can be solved using 
BEM has received considerable mention in the literature during the past decade. Apart from 
the reduced dimensionality and the need for no special domain discretization, other derived 
advantages include: 

1. the ability to handle infinitely large domains: 

2. a much reduced coefficient matrix: 

3. the ease with which singularities are handled: 

4. the restriction of the discretization errors to the boundary, so that the solution is as 
good as the description of the boundary geometry: 

5. the robustness when complex geometries are involved: 

6. the ability to find solutions a posteriori at desired points, not at nodes predetermined 
by the domain discretization: 
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7. the great latitude in solving transient problems by a) using the appropriate time depen- 
dent fundamental solution in the formulation: b) applying the technique in a transform 
domain ( e.g .. Laplace or Fourier): or c) using a time marching procedure. 

Efforts at applying boundary integral techniques to nonlinear problems are quite recent. For 
over a decade the main focus was largely on li near problems such as potential flows (see 
e.g.. Brebbia et al.. [1984]. Liggett A: Liu [1981]). BEM formulations for nonlinear (e.g., 
Lafe et al.. [1981]: Lafe cO Caban [1990]) or those in heterogeneous continua Cheng 

[19&4]: Lafe et oL, [1987-1992] have relied heavil y on iterative methods which still require 
some domain integration. The Dual Reciprocal technique creates a major path for exploiting 
the advantages of BEM to solve nonlinear problems such as convective flows, No domain 
integration is involved when the Dual Reciprocal approach is followed. 

The original credit for Dual Reciprocal BEM concept goes to Savdini k Brebbia [1 98*2] 
who first suggested an innovative approach for transforming domain integrals to the bound- 
ary. However, until recently, prior investigations (see r.r/.. Brebbia [1991]; Partridge et oL w 
[1992]) did not make use of a complete set orglobal functions. A series of local radial func- 
tions were utilized. This made convergence difficult or impossible for a class of nonlinear 
problems. This author and his co-workers (see Cheng ft nl.. [199-1]) have recently derived a 
set of complete coordinate functions which have been tested on a family of strongly nonlinear 
PDEs. Excellent results have been obtained with the complete set. This work opens the 
door to the application of BEM to a wide spectrum of complex flow problems. 

In this paper, we present the full formulation of the Dual Reciprocal Boundary Element 
Method (DRBEM), for incompressible convective flows. 

GOVERNING EQUATIONS 

Let the flow region is represented by Q and t he boundary is T. The pertinent flow equations 
are: 

• Continuity Equation 

V • v = 0 (1) 

• Conservation of Momentum 

ih/ ] I 

— + (v.V)v = --V,,+ -V.7 + g (2) 

where v is the velocity, p is pressure, g is the gravitational accelerator vector, r is the 
viscous stress tensor. If p is viscosity, then for a Newtonian fluid, r is expressible in 
the form: 

T = pVv 


Dimensionless Equations 

Let L = characteristic length scale. ~ = mean velocity, and >/ is the elevation of the point 
(x). We can define the following dimensionless variables: 

x. = x/L, (3) 
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v, = v/r 

p- = (p + pp)/ (n~ 2 ) 

L = fr/L 

With these, the above conservation statements can be made dimensionless: 

Vv. = 0 

^7 + (v -‘ V)v - = - v, '- + 7T vV 


(4) 

( 3 ) 

( 6 ) 


(T) 

(S) 


where 


R f — pvLjp = Reynolds Number (9) 

The governing equations can be rearranged and written in the pseudo- Poisson form: 


= F(x..t.) 


where 


and 


$ = 


v, \ elocity 
p. Pressure 


F = 


li, l<)v,j ()t „ + ( v, • V)v„ + V/),] Velocity Equation 


( 10 ) 
( 11 ) 

( 12 ) 

—V • [(v • V)v] Pressure Equation 

The pressure equation is obtained by introducing the. continuity equation into the divergence 
of the momentum equation. Note that in the velocity equation. <t> and F are vectors with 
two (for 2-D and axi-synunetrie problems) or three (for .'$-D problems) components. We will 
now drop the * prefix in the dimensionless variable's, for convenience. 

For most flow problems the boundary conditions will generally consist of three types: 


• Dirichlet Boundary (F*) 

• Neumann Boundary (F y ) 


<t> = <k 


<)$ 


where dQ/dn — V<!> • n. anil n is the unit vector normal to the boundary. 

• Mixed (r.u) 

(,"($. V$.x. /) = 0 

where (, is some specified function. A free-surface will be an example of the third. In most 
iterative schemes it is usual to recast the M ixeil boundary condition in the form of either 
the Dirichlet or the Seumaun types. 
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BOUNDARY INTEGRAL EQUATIONS 

We will use equation (10) as the representative PDE in developing the integral equations. 
If fictitious sources of strength tr are distributed around T. equal ion (10) can be converted 
into the integral expression (sec Jaswon Synim [1077]): 

$(x) = f ic(x')g(x.x') dx' + f F(x")g{x. x") dx” (13) 

•'r -Iq 

where g is the free-space Green’s function which must satisfy: 

V 2 g{x.x') = fi[x. x') (14) 

where S is the Dirac delta function applied at a point x' and felt at x. The closed form 
solution to equation (11) is {Greenjierg [1971]): 

\nr/2~ in two-dimensions 

g(x.x’) = < ( 15 ) 

1 /( 4 —;• ) in three-dimensions 

in which r — |x — x'|. The last term in equation (13) represents a domain integral. To 
convert this into an integration on the boundary we introduce the Dual Reciprocal concept 
( Cheng et a 1 . . [1993]). 

DUAL RECIPROCAL TECHNIQUE 

Consider nj points on T and in f>. We introduce a family of coordinate functions Mj(x) 
( j = 1. 2, • ■ ■ n j) such that : 

F(x)*£; 1 } M,(x) (16) 

where 3j are coefficients to be determined by rollueat ion. We assume for each function 
Mj{x ). there exists an associated function *I^(x) sitrh that: 

V 2 ^ / (x) = M t (x) (17) 

It can be shown (Cheng A: Oua/ar [ 1 95)3] ) that for a two- dimensional problem for which 

Mj = .r m y n the function tyj is given by: 


q>-l + ^ + l 

‘ im + 2k)l(»-2k+2)\ 


t \ ) (m — 2A’+ > | 4- }! 


where the square brackets in the upper limit of the summation denote the integer part of 
the argument. Solutions for other possible families of coordinate functions are presented in 


Table 1. 
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Table 1 




C 

Dim. 

.17, 

*7 

V 2 

2D 

cos(/?.r) cos (/»//) 

zlh 

( ;t 2 + m 2 ) 

V 2 

Axi 

A’o(jtr) cos (kz) 

-y 

(n J + le ) 

V 2 

3D 

cos (n.r) cos(/n//) cos( kz) 

( + 2 + k 2 ) 

V 2 

2D 

^ (njr+m.v) 


V 2 

Axi 

K 0 (nr)e kz 

Ih 

{n 2 + k 2 ) 

V 2 

3D 

f (nx+ni.v+k?) 

1L 

V 2 - A 2 

2D 

cos(n.r) cos (my) 

-Vb 

( n - + m 2 — A 2 ) 

V 2 - A 2 

Axi 

I\o(nr) cos(k-z) 

- v b 

(r ,-'+*.2-.V) 

V 2 - A 2 

3D 

cos(n.r) cos( my) cos(7 - c) 

~'b 

(„2 + „,-> + fc2_\2) 

V 2 - A 2 

2D 

t (»r+m.v) 

Hi 

( + — \ 2 } 

V 2 - A 2 

Axi 

Ko(»r)r ts 

1L 

(„2 + A.2_.\7) 

<1 

! 

3D 


1/, 


( n 2 +m 2 +k 2 — \ 2 ) 


in which Ao is the zeroth order modified Bessel function of the first kind. When equations 
(16) and (17) are used in (10). and we distribute the fictitious sources on T we can obtain 
the "pure' boundary integral equation: 

ir{x’)g(x.x') dx' + *,-¥,-(x) (19) 

i= i 

An expression for the gradient of <t>. which is required in equation (12) can be obtained from 
equation (19) as: 

E ( 2 °) 

7 = 1 


V$(x) = f u-(x')Vg{x. x) dx' + 
r 
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( 21 ) 


The normal derivative d$/On = V$ • n is given by: 



uix')^(x.x') <lx' + : ^(X) 


Assuming 3j ( j = 1. 2 * * * nj) are known, the only unknown in equations ( 19) and/or (21) 
is the source strength distribution ir on T. The trick is to start with a trial distribution of 
i r (x) and to obtain the coefficients 3j ( j = 1.2- — nr) by collocation using equation (16). 
When applied to all nj selected points the result of the collocation is the matrix equation: 


n T 


H-UM = Ft i = 1-2. — /?x 


( 22 ) 


where A/y — Mj(x,) and F, = F(x, ). The matrix system (22) is also expressible in the form: 


MJ = F 


(23) 


which can be inverted to give: 

3 = M -1 F (24) 

Once 3 has been determined, equations (19) and/or (21) are then combined with the 
prescribed boundary conditions and solved for ic on Y. A better estimate of F is then 
obtained by using equations (19) and (20) in (12). The solution process continues until a 
specified convergence criterion is satisfied. 


DISCRETIZATION 

We subdivide the boundary into tu, elements. Let A*.(x) (k = 1.2. •••»)$) represent the 
boundary shape functions describing the distribution of tr on I". By selecting each of the n j, 
boundary points as successive origins oi integral ion. equations ( 19) and (21 ) can be assembled 
into the svstem: 


"b 

YY a ile w k = l>,. i = 1 . 2 . 


where 


A = l 


«ik = 


bi = 


"b 

x, e r$ 
x, 6 1 q 

x, s r* 


f\\ A "*(x / )$r(x'.x,) 3x' 

Sr, A/ (x' )dgf On ( x'. x, ) fix' 

AfxA - V nr LU/ 

M/Mx;)- Z"Lx ijM-jOn x, € Y q 

tih). Symbolically equation 


(25) 

(26) 
(27) 


Therefore, we have ni, equations to determine u\. (k =1.2. 
(25) can be written in the alternative form: 


aw = b (28) 

which can be inverted to give: 

w = a~‘b (29) 

The whole process boils down to the iterative solution of equations (24) and (29). with 
repeated updating o( F using (12). At each time level /. the iterative steps are: 

1. Start with a trial F (/.r.. F, values for i^.= l. 2. • • • nr). 
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2. Obtain from equation (24). 

3. Obtain w using equation (29). 

4. Ise discretized forms of equations (19) and (20) to compute at all nj points. 

This provides a better estimate for F. 

5. Go back to Step 2 il convergence condition is still unsatisfied. 

Note that the matrix inversions in equations (24) and (29) need only be performed once, 
for fixed boundary problems. The vectors w and 4 are the quantities whose values change 
during the iterative process. Once convergence is reached, equations (19) and (20) can be 
used routinely to obtain = (v„.p, ) or the gradient at any point (x) of interest. 


Treatment of Time Derivatives 

We now need to address the treatment of the local acceleration term dv/dt as occurs when 
equation (11) is written for the velocity. We discuss two efficient approaches for handling 
this term. The first is based on the use of a time- dependent fundamental solution. The 
second utilizes a simple time- marching procedure. 

Time-dependent Fundamental Solution 

Equation (11), written for the velocity {/.(.. $ = v). can be re-arranged into the alternative 
form 

£$(x./) = F(x.l) (30) 

where 


c = 

of 

F = R. [(v • V)v + Vp] 
* 

The boundary integral equation in this case is 



where the functions g and 'P, must respectively satisfy the following PDEs 


Cgix.t.x'.t') = (‘'(x. /. x\ /') 
£tf(x./) = Mj(x.t) 

The closed-form solution for (32) is (see Greenberg [1971]) 


g{x.t.x.l ) = { 


| * ex P ip-c) } in two-dimensions 

1 ^I'T/ 1 ex P (— } in three-dimensions 


(31) 


(32) 

(33) 


(34) 
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where H is the Heaviside function. 

The extra computational effort required here is the time integration. / 0 f (-)< at each time 
level. However, this approach has the major advantage that no difference approximation is 
introduced in the evaluation of the time derivative, and the exact time-dependent funda- 
mental solution is utilized in the integral equation. 


Time-marching Procedure 

In this approach we assume the time derivative can be approximated by the difference equa- 
tion 

dv v — Vu 
~dt ~ It 

where Vo is the velocity at a previous time level. The velocity equation (30) is still valid but 
the differential operator C and forcing function F now become 


C 

F 


v 2 (-) — -vt(') 


ft. 


A / 

+ ( y 0 ■ V)v 0 -I- V/» 


The boundary integral equation in this case is 


u*(x'. /')</(x. x") fix' + ^2 ij'Vjix.i) (35) 

j=I 

where the functions <j and ty, must respectively satisfy 

Cg{x.x') = t){x.x') 

CV(x.f) = Mj(x.l) 

It is easily shown that 

{ A it in two-dimensions 

^exp{ — in three-dimensions 

The time-marching approach has the obvious advantage of not requiring an explicit time 
integration, as in the first method. Furthermore, for fixed boundary problems, the free- 
space function g need not be calculated at each time level. The iterative scheme can now be 
replaced by the time-marching process. However, the presence of A / and R, in the argument 
of the Green s function creates the need for extreme caution in the numerical evaluation of 
g. As the time step is reduced to improve accuracy, or as the flow moves away from the 
laminar regime, the numerical value of g reduces very rapidly. 



(36) 

(37) 
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CONCLUSIONS 

New concepts in boundary element modeling provide excellent approaches for solving con- 
vective transport problems. A formal determination of the advantages inherent in the new 
DRBEM formulation can easily be accomplished. In general, we expect optimum nj (=total 
number of collocation points selected to evaluate coefficients 3) to be of the same order as 
the number, nb, of boundary elements. Therefore, the largest matrix size in BEM will be of 
order x rib = By comparison, the domain techniques, because of the need to discretize 
the entire region, will produce a global matrix size of order nj x = nj- The advantage 
in terms of storage requirements is obvious. Moreover, the much reduced size of the global 
matrix has a more pronounced advantage in total ('PI' time. 

For example, in a 2-D flow problem discretized with 100 boundary elements, the matrix 
size using the domain methods (if no special consideration is given to matrix bandedness) 
will be 10 4 larger than that of BENL Even with the sparseness of the global matrix taken into 
account in the domain methods, a boundary element approach still enjoys a size advantage 
proportional to the bandwidth of the matrix. The computational advantage in 3-D is more 
significant because of the much increased number of meshes in the domain techniques. 

These computational advantages are key to effective modeling of convective flows. The 
compactness of BEM matrices allows for a greater freedom to experiment , even on computers 
of average memory! DRBEM provides an excellent platform for optimizing system geometry. 
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